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Abstract 

The vortex glass transition in the presence of columnar defects is studied by 
Monte Carlo simulations of a vortex loop model, suggested by the analogy 
to the T = superconductor-insulator transition for dirty bosons in (2+l)D. 
From finite-size scaling analysis of the I-V characteristic we find two dynami- 
cal exponents describing the non-equilibrium behavior. We obtain z± = 6±0.5 
and Z|| = 4 lb 0.5 when the current is applied perpendicular and parallel to 
the columnar defects respectively. 

PACS numbers: 74.60.Ge, 67.40.Fd, 74.70.Vy, 74.70.Mq 
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Strong thermal fluctuations and quenched disorder give properties to high temperature 
superconductors not observed in conventional type II superconductors. The irreversibility 
line separates two regions in the mixed state with distinct properties of the penetrating 
magnetic field lines. Above the irreversibility line, vortex lines are in a liquid state due 
to the strong thermal fluctuations. Matthew Fisher [|T| has suggested that in the presence 
of quenched disorder there is a phase transition into a new superconducting state at the 
irreversibility line, called the vortex glass One signature of a vortex glass transition 
is universal scaling of the nonlinear current- volt age (I-V) characteristic This has been 
verified in several experiments 0-0], and in Monte Carlo simulations |]8HTl • 



The original vortex glass model assumes pointlike impurities. In a recent experiment 



Civale et al. |T2| found a remarkable shift upwards in the irreversibility line of several Tesla 
upon irradiating a sample of YBa2Cu307 by heavy ions. The damage tracks from the ions 
form permanent "columnar defects" with linear dimension comparable to the sample size. 
These columnar defects strongly enhance flux pinning in the superconductor. 

Nelson and Vinokur |T^ derived scaling relations for the I-V characteristic for the case 



of columnar defects, using an analogy to the Bose glass []I1|,|T3[. In this picture, the glass 
transition is analogous to the zero temperature superconductor-insulator transition of dirty 
bosons in two dimensions. The magnetic flux lines through the 3D bulk superconductor 
correspond to world lines of bosons in the (2+l)D path integral formula for the partition 
function. Static disorder for the (2+l)D boson system corresponds to columnar defects for 
the 3D superconductor. The superconducting phase corresponds to the insulating localized 
Bose glass phase for dirty bosons, and the resistive vortex fluid phase corresponds to the 
superconducting phase of the dirty bosons. 

A question of principal interest is whether a stable glass phase can actually exist in 3D. 
This issue was studied by Monte Carlo (MC) simulations of gauge glass models with point 
disorder P,p|,p!6| . These found that the lower critical dimension is somewhat smaller than 



c? = 3, but so close to = 3 that it is hard to resolve unambiguously. On the other hand, the 
existence of the glass phase in the presence of columnar defects is clear due to the analogy 



with the dirty boson problem which has a phase transition which is now well established 
T^,[18|. Furthermore, in these calculations screening of the flux lines, and the anisotropy 



due to the preferred direction picked out by the applied magnetic field, were not considered. 

In this Letter we report Monte Carlo (MC) results for the nonlinear I-V characteristic at 
the glass transition of a model 3D high temperature superconductor with columnar defects. 
Screening of the flux lines is taken into account by neglecting long-range vortex interactions 
altogether and assuming only a short-range repulsion. Our model effectively assumes that 
the magnetic field fluctuations diffuse more rapidly than the vortices. The more difficult 



case of long-range forces and columnar defects is being considered elsewhere ||T9 |. 

Let us briefly discuss scaling of the nonlinear I-V characteristics at the glass transition 
||l|J13[| . Any configuration of wiggling vortex lines may be described as a constant background 



of straight lines plus closed loops. The characteristic size of the loops defines the correlation 
length ^, which diverges as ^ ~ \T — Tg\~'^. The correlation time is assumed to obey standard 
dynamical scaling: r ~ where z is the dynamical exponent which will be computed 
below. Columnar defects make the system strongly anisotropic, and this is reflected in the 
divergence of the correlation length. From the analogy to the dirty boson problem we infer 
that the correlation length in the direction of the columns scale as ~ where zq = 2 
is the quantum "dynamical exponent" of the dirty bosons at T = in 2D ||T8| . 

Consider voltage fluctuations across a correlation volume of size x ^_|_ x ~ ^'^^'^ . 
The linear response electric field E, the current density /, and the resistivity p are related by 
E = pi. From the Josephson and Nyquist relations the resistance R across the correlation 
volume scales as i? ~ In the case of a transverse current applied in the direction 

perpendicular to the columnar defects, the linear resistivity becomes p = R$,±$,\\/$,± ~ ^ , 
and E ~ ^^"^ ^I. We also consider the case of an applied longitudinal current. Here we find 
p = R^l/^^^ ^ el"''''', and E ~ ^I'^^^-'I. 

We now have to find the nonlinear I-V characteristics. At finite current density /, vortex 
lines fall out of equilibrium due to the Lorentz force. The work against the Lorentz force to 
create a closed loop costs energy E ~ ±AI, where ± is the orientation and A the area of 



the loop. For large enough loops the energy shift reaches ksT. This defines a characteristic 
"current length" C,i beyond which correlations are destroyed In an isotropic system we 
would define ^fJ = 1 (dropping a factor of ksT). With columnar defects and a transverse 
current we instead take {^±^\\)il ~ ^j^^'^I = 1. The nonlinear I-V characteristic is now 
expected to have the form V ~ with exponent y± = {1 + z)/{1 + zq), and, in general, the 
electric field is expected to scale like E ~ F±{I^^^'^), where and F_ are universal 

scahng functions valid for T >Tg and T < Tg, respectively [^]. In the case of a longitudinal 



current we instead have {^jjil ~ ^jl = 1, and following the same derivation gives V ~ /^n 
with 2/11 = {ZQ + z)/2, and E ~ C^''''''^ F^^i^e^)- 

For a finite system at T = Tg the diverging correlation length is cut off by the system 
size L. With transverse / we therefore expect E ~ Thus a plot of EL^^^ 

versus IL^^^'^ is expected to give a universal, system size independent curve. Notice that 
the meaning of the scaled quantities is very natural. We have EL^^^ ~ V/{V'^Ye^^ where 
V = EL, (l^^)e^^ is the rms equilibrium voltage fiuctuation, and IL^'^^'^ ~ {L/^iY^^'^ . 

For our simulation we want a minimal model that describes the collective behavior of 
interacting vortex lines in the presence of thermal fiuctuations and columnar defects. We 



choose the dirty boson action 



(iH = (3Y.{b\^)-v{v^)Uv)], (1) 

where [3 = l/fc^T, J = (J^:, Jyi Jz) are integer variables defined on the links of the lattice, 
taking values from — oo to oo. These integer variables represent the vortex lines, for example, 
J2:(r) = ±1 means that one vortex line with orientation ±1 goes from r to r + x. The vortex 
"current" is constrained to be divergenceless. The first term in H acts as a string tension 
and as a short-range repulsion. The second term in H describes the columnar defects, f (r^) 
is a random site energy with uniform probability distribution on the interval [0, 1]. 

Due to the anisotropy from the columnar defects, the correlation length in the direction 
of the columns diverges at Tg as ~ with zq = 2. To enable finite-size scaling we must 
use lattice sizes that scale correspondingly. We use lattices of size L x L x with periodic 



boundary conditions in all directions. We use an applied magnetic field such that there is a 
fixed density of / = 1/2 of vortex lines in the direction of the columnar defects, corresponding 
to half filling in the dirty boson problem. The transition point Tg = 0.248 ± 0.002 for this 
model is already known from our previous work on the superconductor-insulator transition 



Our MC algorithm consists of attempts to insert closed loops of vortex lines with random 
sign on the plaquettes of the lattice. The algorithm is ergodic in the sense that all possible 
configurations of closed loops are accessible. Non-equilibrium effects due to the Lorentz 
force from the finite current are included by biasing the acceptance of plaquette moves. The 
energy change of a move includes a term of the form Ae = (V x J) ■ I, where the discrete 
curl gives the orientation of the plaquette loop. 

We do a thermal average over vortex loops for each realization of the columnar disorder, 
followed by a disorder average. As an estimate of the required number of MC update sweeps 
in the thermal averages we use the time to return to steady state after a reversal of the 
current. The voltage is measured over up to 10^ sweeps, with up to 10^ initial sweeps 
discarded to reach steady state. For large currents a smaller number of sweeps is sufficient. 
Also the number of disorder realizations necessary to get small statistical error depends on 
the size of the applied current. At large enough current, so that ^/ <^ L, the system will 
self-average and a small number of disorder realizations is enough. At small I we use up to 
20 realizations of the disorder. 

The voltage is given by the rate of phase slip across the system. This is measured by 
keeping track of the number of plaquette moves in the plane perpendicular to the applied 
current. The electric field is proportional to 



where Q = LxLyLz is the volume of the system, N± is the number of accepted plaquette 
moves with positive/negative orientation in the plane perpendicular to the applied current, 
and the time derivative stands for the change per sweep. The bracket ( ) denotes thermal 



m 
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average, and [ ] denotes quenched disorder average. Equation (||) implicitly assumes heavily 
overdamped dynamics so that MC time can be equated with real time. The update order 
in a sweep through the lattice is random rather than sequential, in order to avoid possible 
systematic error at high current where the acceptance rates are high. 

Monte Carlo results for the nonlinear I-V characteristic at the glass transition are shown 
in Fig. |l|. In (a) the current density / is transverse to the columns, and in (b) / is longitudinal 
20| . The dynamical exponent z was adjusted until all data for different lattice sizes collapse 



on the same curve. Surprisingly, this happens for different values of z for a transverse and 
a longitudinal current. In (a) z has been set to z± = 6, and in (b) z has been set to 
z\\ = 4. The MC data has three different regimes: (1) For small current, where ^/ > L, 
the MC data points approach the dashed line, whose slope is 1. This is a finite size effect: 
when the current length is cut off by the system size L the response becomes ohmic. (2) 
In the nonlinear I-V scaling regime, 1 < ^/ < L, the MC data nearly coincides with the 
solid line, whose slope is given by the scaling prediction for the exponent y in V ^ P. 
In (a) = {1 + 2;_l)/(1 + zq) = 7/3, and in (b) y\\ = {zq + z\\)/2 = 3. (3) At higher 
current, with ^/ < 1, scaling breaks down as expected (data not shown). Error bars on the 
MC points represent the standard deviation of the voltage fluctuations between different 
disorder realizations. 

What is the significance of the scaling plots in Fig. [l|? The critical temperature Tg was 
determined independently 0], and is not an adjustable parameter. The figure shows that 



by adjusting only a single parameter, the dynamical exponent z, we simultaneously get 
two things: (1) all MC data collapse on the same universal curve, and (2) the curve is a 
remarkably good power law over six decades with an exponent which is consistent with the 
predictions of the scaling analysis. For other choices of z both these properties quickly break 
down. From the rate at which the data collapse breaks down, we estimate = 6 ± 0.5 and 
z\\ = 4 ± 0.5. One would have expected, at least in equilibrium, that a single dynamical 
exponent would suffice, but our simulation shows that in this non-equilibrium case we need 



two Attempts to use a single exponent z 5 give a distinctly poorer fit to the data. 



As a test of universality we changed the aspect ratio of the lattices from = to both 
and L^. We also varied the applied magnetic field from / = l/2to/ = l/4. These did 
not change the values of the dynamical exponents. We have not yet tested lattice structures 
other than simple cubic. 

The following simple argument suggests how the two different dynamical exponents z± = 
6 and z\\ = 4 might arise. First consider an isotropic system, and study the motion of a 
blob of vortex lines of diameter ^. As a crude approximation we view the blob as a polymer 
containing ~ segments. The time it takes each segment to random walk a distance ^ 
scales as ti ~ C,'^. The diffusion coefficient for the center of mass (CM) is down by a factor 
of A^. Hence the time for the CM to move a distance ^ scales as t ~ and thus z = 4. 
This argument neglects correlations in the motion of the segments of the vortex lines. The 
result coincides with the mean field result . When a finite current is applied, the role of 



the diverging correlation length ^ is taken over by the finite current length which in the 
isotropic case is defined by ^j'^^I = 1. The time for the blob to move is now t ~ 
so again z = 2(1 + 1) = 4. In the anisotropic case with columnar defects and a transverse 
current, the current length is defined by ^j'^^'^ I = 1. This is now the characteristic length 
that the blob has to move. If we assume that the time for the CM of the blob to move 
is now t ~ ^^(^+^<5)^ = 2(1 + zq) = 6. With a longitudinal current we instead 

have = 1, and the correlation time is t ~ so z\\ = 2(1 + 1) = 4. If our 

assumption is justified, this argument could explain why the non-equilibrium correlation 
time is different when the current is applied in different directions. Exactly what happens 
when the applied current goes to zero is not clear from this argument, however one normally 
would expect a single relaxation time. We note in passing that a measurement of the voltage 
noise spectrum at the critical point should be able to independently determine the predicted 
current dependence of the characteristic frequency scales in the system. 

How do our results compare with experiments on systems without deliberately intro- 
duced columnar defects? Experiments PBI on single crystals of YBCO give z = 4.5 ± 1.5 



and u = 2.0 ± 1, by fitting data to point-disorder vortex glass scaling laws. In these experi- 



ments correlated disorder might be present in the form of dislocations, twin boundaries, etc. 



If, following Nelson and Vinokur |T3[, we instead fit to the Bose glass scaling laws used in 
this paper for columnar defects, we find z = 7 ± 2 and u = 1.3 ± 0.5. These are roughly 
consistent with our values z± = 6 ± 0.5 and u = 1.0 ± 0.1, where u was calculated in pre- 
vious equilibrium simulations of dirty bosons . Simulations of point and twin-boundary 
disorder is underway. 

In summary, we report Monte Carlo simulations of the nonlinear I-V characteristic of 
a model high-temperature superconductor with columnar defects. We find two dynamical 
exponents describing the non-equilibrium behavior, z± = 6 ±0.5 and z\\ = 4 ±0.5 in the case 
of an applied transverse and longitudinal current respectively. Experimental measurements 
in systems with deliberately introduced columnar defects would be highly desirable. 

Many helpful discussions with D. Nelson, M. P. A. Fisher, D. Stroud, D. Huse, A. Sudb0, 
A. P. Young, N. Goldenfeld, and A. Dorsey are gratefully acknowledged. MW is supported 
by NFR grant F-DP 9902-303. SMG is supported by DOE grant DE-FG02-90ER45427 and 
NCSA DMR-910014N. 
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FIGURES 

FIG. 1. Log-log plot of Monte Carlo results for the nonlinear I-V characteristic at the vortex 

glass transition with cohimnar defects. Plotted in (a) is MC data for EL^^^ versus IL^'^^'^, where 
E is the electric field, and the applied current density / is perpendicular to the columnar defects. 
Plotted in (b) is EL^Q'^^ versus /L^, where the current density / is parallel to the columns. The 
dynamical exponent z has been adjusted until all data for different system sizes collapse on the 
same curve. In (a) z has been set to z± = 6, and in (b) z\^ = 4. For other values of z scaling 
quickly breaks down. The solid straight line in (a) has slope y± = {1 + z±) / {1 + zq) = 7 /3 given 
by scaling (see text). The solid straight line in (b) has slope y\\ = {zq + z\\)/2 = 3. The dashed 
straight lines in (a) and (b) both have slope 1, and indicate the ohmic region at small current. 
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